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ABSTRACT 


Two  Fortran  subroutines  (BICUB1  and  BICUB2)  which  perforin 
bicubic  spline  interpolation  of  a  tabulated  function  of  two 
variables  are  described.  Given  the  values  X  (1) , . . .  ,X  (N)  and 
Y(l) ,. . . ,Y(M)  of  two  independent  variables  and  the  corresponding 
function  values  |U(I,J)-f  (X(I) ,  Y(J))|,  I-1,...,N  and  J«1,...,M 
and  certain  normal  derivatives  (optional)  along  the  boundaries 
of  the  x-y  mesh,  BICUB1  estimates  the  derivatives  f  ,  f  ,  and 

X  7 

f  at  each  (I,J)  mesh  point.  If  the  normal  derivatives  along 
xy 

the  mesh  boundaries  are  unknown,  BICUB1  estimates  them  using  a 
moving  third  order  two  dimensional  Lagrange  interpolating 
polynomial.  Given  the  coordinates  (XPT.YPT)  and  the  derivatives 
calculated  by  BICUB1,  BICUB2  obtains  the  coefficents  of  the 
bicubic  polynomial  for  the  rectangular  region  of  the  mesh  con¬ 
taining  (XPT,YPT)  and  estimates  the  functional  value 
UPT«f (XPT,YPT) .  In  effect,  the  routines  pass  a  twice  continuously 
differentiable  piecewise  bicubic  polynomial,  u(x,y)CC  ,  through 
the  given  functional  values. 


PROBLEM  STATUS 


This  is  a  final  report  on  one  phase  of  a  continuing 
problem. 


AUTHORIZATION 

NRL  Problem  S01-38 


ii 


1.0  IDENTIFICATION 


1.1  Title 

Bicubic  Spline  Interpolation 

1 .2  Identification  Name 
El -NRL-BICUBIC 

1.3  Classification  Code 

El 'Interpolation  and  Approximations,  Curve 
Fitting 

1.4  RCC  Identification  Number 
E1OO10OO 

1 . 5  Entry  Points 

BICUB1 

BICUB2 

1.6  Programming  Language 
Language:  3600/3800  FORTRAN 
Routine  type:  Subroutines 
Operating  System:  DRUM  SCOPE  2.1 

1.7  Computer  and  Configuration 
CDC  3800 

1 .8  Contributor  or  Programmer 
John  J.  Cornyn* 

Information  Processing  Systems  Branch  (Code  5493) 
Communications  Sciences  Division 


formerly  with  the  Large  Aperture  Systems  Branch,  Code  8160, 
Acoustics  Division 
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1  * 


1.9  Contributing  Organization 


NRL  -  Naval  Research  Laboratory, 

Washington,  I).  C.,  20375 

1.10  Program  Availability 

1.10.1  Submittal :  Program  write-up,  Fortran 

source  deck,  source  listing 

1.10.2  On  File:  RCC  Program  Library 

1.11  Verification 


Several  third  degree  polynomials  were  used  to 
test  BICUBIC;  answers  were  good  to  at  least  nine 
significant  figures.  Higher  degree  polynomials 
were  also  used.  Then,  as  expected,  the  results 
did  not  compare  as  well  with  the  true  values. 

1.12  Date 


26  February  1973 
2.0  PURPOSE 

2.1  Description  of  Routines 

Let  the  values  of  a  function  u(x,y)  over  a  two 

dimensional  domain  be  given  at  the  mesh  noints 
(xi,yj)  where  i  =  1 ,  ...,N;  j  «  1,  . ..,  M. 

(1)  The  first  problem  considered  is  the  estimation 
of  the  normal  derivatives  along  the  boundaries 
of  the  mesh  assuming  they  are  unkr  >wn.  In 
Figure  1,  squares  designate  locations  at  which 
one  needs  to  know  the  normal  derivatives  with 
respect  to  x,  p..  =  u  (x.,y.).  Circles 

1 J  A  1  j 

designate  locations  at  which  one  needs  to  know 
the  normal  derivatives  with  respect  to 
y,  =  Uy(x^,  Vj).  Squares  imbedded  in 

circles  designate  locations  where  the  normal 

derivatives  with  respect  to  both  x  and  y, 

S..  =  u  (x. ,  y.),  are  required,  in  addition 
xy  l  j 

to  p^  and  q^j .  A  solution  to  this  problem 

will  be  given  in  the  form  of  subroutines 
EDGES  and  LAGRAN. 
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(2)  The  second  problem  is,  assuming  the  normal 
derivatives  along  the  boundaries  as  discussed 
above  have  been  given  or  estimated,  to  fit  a 

2 

"smooth  function"  u(x,y)C  C  (twice  continuously 
differentiable)  through  these  values. 

The  bicubic  spline  interpolation  routines 
BICUB1 ,  BICUB2 ,  GETBP,  AND  SOLVTT 
described  later  implement  a  bicubic  spline 

interpolation  technique ^  [l]  (see  Section  3. IS) 

which  yields  a  piecewise  bicubic  polynomial 
u(x,y).  This  function  is  defined  in  each 
rectangular  cell, 

Rij  :  xi-l-x-xi  5  yi-l-yiyi  »  0) 

of  the  grid  as 

3 

u(x,y)  -  Cj.(x,y)  -£  T*J(x-xi.1)*  (y-y).1)n  (2) 

m,n»0 

where  (x.yjcRjj. 

Individual  Subroutine  Functions 

BICUB1  -  A  subroutine  for  calculating  the 

normal  derivatives  at  each  mesh  noint. 

BICUB2  -  A  subroutine  for  interpolating  a  value, 
u(x,y),at  any  point  (x,y)  within  the 
region  subtended  by  the  mesh. 

EDGES  -  A  subroutine  for  estimating  the  re¬ 
quired  normal  derivatives  along  the 
boundaries  assuming  they  have  not 
been  given,  using  a  moving  third 
order  two  dimensional  Lagrange 
interpolating  polynomial. 


Note  that  Eq.(10)of  [l]  contains  a  typographical  error  and 
should  read 

T 

A ( Ax i _ i ) ^i j  A  (Ayj_i)  °rij»  where  T  indicates  a 

transpose.  Also,  the’ 8th  character  on  line  4,  j\2 1 2  of  I’ll 
should  be  I  and  not  1.  L J 
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GETBP 


■  A  subroutine  for  calculating  the  two 
dimensional  difference  arravs  b  and  b* 
of  Eq. CIS)  of 


LAGRAN  -  A  subroutine  for  determining  the  value 
of  a  two  dimensional  Lagrange  inter¬ 
polating  polynomial  of  degree 
m  in  x  and  n  in  y  and  its  derivatives 
with  respect  to  x,  with  respect  to  y, 
and  with  respect  to  both  x  and  y  at 
any  intersection  of  a  two  dimensional 
mesh  defined  by  Cro+1)  levels  of 
x^,  i  *  a,  a+m  and  (n+1)  levels  of 

Yy  j  s  c»  •••»  c+n.  See  Eq.(3). 

SOLVIT  -  A  subroutine  for  solving  a  linear 

system  using  Gaussian  elimination  as 
illustrated  in  Eqs.(15)  and  (16)  of 

[i]  • 

2.2  Problem  Background 

In  the  development  of  models  of  certain  physical 
phenomena  it  is  frequently  useful  to  obtain  a 
smooth  functional  representation  of  a  quantity 
which  is  known  at  only  a  discrete  set  of  points 
over  a  two-dimensional  domain.  Frequently,  it  is 
required  that  this  function  have  continuous  first 
and  second  derivatives.  Once  such  a  representation 
is  obtained,  it  is  possible  to  differentiate  and 
integrate  it  in  closed  form. 

A  major  problem  with  a  Lagrange  interpolating 
polynomial  defined  over  a  N  x  M  mesh  is  that  it 
must  maintain  (N-2)  in  x  and  (M-2)  in  v  continuous 
derivatives  and  still  pass  through  all  of  the  data 
points.  These  requirements  can  and  generally  do 
lead  to  many  large  and  unrealistic  mountains  and 
valleys  in  the  interpolation  surface,  i.e.,  large 
interpolation  errors.  This  problem  is  significantly 
reduced,  but  not  completely  eliminated  by  the 
bicubic  spline  (See  Section  3.13). 


*  r 
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A  secondary,  but  still  significant,  problem  is 
that  when  a  large  number  of  points  is  used,  the 
evaluation  and  calculation  of  the  Lagrange 
interpolating  polynomial  is  costly  and  unreliable 
([6]  ,  p.  231).  For  further  discussion  see  £8]. 


3.0  USAGE 


3.1  Calling  Sequence  or  Operational  Procedure 

BICUB1  (N,M,NDFAULT,X ,Y,U,P,Q,S,MAX  ,W1  ,W2  ,W3  ,W4 
W5,W6,W7) 

BICUB2(XPT,YPT,UPT,NERR0R,N,N,X,Y,U,P,Q,S) 

3 . 2  Arguments,  Parameters,  and/or  Initial  Conditions 

N  --  number  of  x  points  at  which  the  function 
was  observed. 

(N  i  4)  .  TYPE  INTEGER. 

M  --  number  of  y  points  at  which  the  function 
was  observed. 

(M^4).  TYPE  INTEGER. 

NDFAULT  --a  parameter  that  must  be  set  to  1  if 

subroutine  BICUB1  is  to  call  subroutine 
EDGES  to  calculate  the  "required"  normal 
derivatives  along  the  boundaries  of  the 
mesh.  If  NDFAULT  is  not  set  to  1,  subroutine 
BICUB1  assumes  the  normal  derivatives  for 
the  boundaries  have  already  been  entered 
into  arrays,  P,Q,  and  S,  by  the  user's 
calling  program.  TYPE  INTEGER. 


The 


"required" 
P(I,J)  for 
Q(J»I)  for 
S(J,I)  for 


normal  derivatives  are: 
1=1  and  N;  J=1  to  M. 

J=1  and  M;  1=1  to  N. 

1=1  and  N;  J=1  and  M. 
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X  --  the  vector  of  distinct  values  of  the  first 
independent  variable  arranged  in  ascending 
order;  (x^  for  i  *1  to  N) .  The  minimum 

length  of  X  is  4  and  the  maximum  length  is 
determined  by  the  amount  of  core  available. 
DIMENSION  X(N).  TYPE  REAL. 

Y  --  the  vector  of  distinct  values  of  the  second 
independent  variable  arranged  in  ascending 
order,  (y^  for  j=l  to  M) .  The  minimum 

length  of  Y  is  4  and  the  maximum  length  is 
determined  by  the  amount  of  core  available. 
DIMENSION  Y(M).  TYPE  PEAL. 

U  --  the  matrix  of  function  values  corresponding 

to  X  and  Y,  i.e.,  U(I,J)  is  u.  .  of  Section  2.1. 
DIMENSION  U(N,M)  .  TYPE  REAL.1-1 

P  --  the  matrix  of  normal  derivatives  with  resnect 
to  x  corresponding  to  X  and  Y;  i.e.,  P(I,J) 

is  U  (X .  ,Y . )  .  If  NDFAULT  is  not  1,  P(I,.T)  for 
g  ^  ^  J 

C=1  and  N;  J=1  to  M)  are  required  input.  If 
NDFAULT  is  set  to  1,  no  values  of  P  are 
required  as  they  will  be  calculated  bv  sub¬ 
routine  EDGES.  DIMENSION  .  TYPE  PEAL. 

Q  --  the  matrix  of  normal  derivatives  with  respect 
to  y  corresponding  to  X  and  Y;  i.e.,  0(.T,I) 
is  Uy(X.,  Yj )  .  Note  thp  inversion  of  the 

J  and  I  indices  in  Q.  Tf  NDFAULT  is  not  1, 

Q ( J , I )  for  (J=l  and  M ; X *1  to  N)  are  reauired 
input.  If  NDFAULT  is  set  to  1,  no  values  of 
Q  are  required  as  they  will  be  calculated  by 
subroutine  EDGES.  DIMENSION  Q(M,N) .  TYPE  PEAL. 


S  --  the  matrix  of  normal  derivatives  with  respect 
to  both  x  and  y  corresponding  to  X  and  Y: 
i.e.,  S(J,I)  is  U  (X.,  Y.).  Note  the 

inversion  of  indices  I  and  J  in  S.  If  NDFAULT  is 
not  1,  s(.J,I)  for  (1  =  1  and  N;  .1=1  and  '0  are 
required  input.  If  NDFAULT  is  set  to  1,  no 
values  of  S  are  required  as  they  will  be 
calculated  by  subroutine  EDGES. 

DIMENSION  S(M,N).  TYPE  REAL. 


fe.. 
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MAX  --  the  greater  of  N  and  M.  TYPE  INTEGER. 


W1 


,W7 


XPT 

YPT 

UPT 


seven  arrays  which  are  used  as  working  areas 
by  subroutine  BICUB1.  Each  of  these  arrays 
must  be  dimensioned  to  MAX  words  in  the  user's 
program.  The  user  does  not  assign  values  to 
these  arrays.  TYPE  REAL. 

the  X  coordinate  of  the  point  to  be  inter¬ 
polated.  TYPE  REAL. 

the  Y  coordinate  of  the  point  to  be  inter¬ 
polated.  TYPE  REAL. 

the  interpolated  value  to  be  obtained.  TYPE 
REAL. 


NERROR  --  an  error  indicator.  If  the  point  (XPT, YPT) 
does  not  lie  within  the  mesh,  NERROR  will  be 
set  to  1  by  BICUB2;  otherwise  it  will  remain 
set  to  0.  IF  NERROR  is  returned  as  1,  the 
interpolated  value  is  set  to  -0.0  and  an 
error  message  is  printed.  See  Section  3.5. 
TYPE  INTEGER. 


3.3  Space  Required  (Decimal  and  Octal) 

3.3.1.  Unique  Storage  (exclusive  of  system  library) 


Subroutine 

Decimal 

Octal 

BICUB1 

1083 

2073 

BICUB2 

570 

1072 

EDGES 

474 

732 

GETBP 

183 

267 

LAG RAN 

588 

1114 

SOLVIT 

129 

201 

TOTAL  3027 

5723 

To  interpolate  a  function  over  an  N  X  M  mesh  the 
user's  program  is  required  to  dimension  the 
arrays  X,Y,U,P,Q,S,  and  W1  through  W7 .  These 
arrays  require  a  total  of 

4NM  +  N  +  M  +  7K  words, 
where  K  is  the  greater  of  N  and  M. 
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3.3.2. 


Common  Blocks 


None . 

3.3.3.  Temporary  Storage 

Once  subroutine  BICUB1  has  been  called  and  the 
elements  of  the  U,P,Q,S  arrays  have  been 
determined,  the  space  occupied  by  the  working 
arrays  W1,W2,...,W7  can  be  used  for  other 
purposes.  That  is,  the  7K  term  of  the  above 
equation  may  be  deleted, 

3.4  Messages  and  Instructions  to  the  Operator 
None. 

3 . 5  Error  Returns,  Messages,  and  Codes 

Subroutine  BICUB1  may  print  out  the  following  error 
messages : 

"ERROR  -  THE  Y  VECTOR  IS  NOT  ARRANGED  PROPERLY. 

ERROR  DETECTED  BY  BICUB1 . 

THE  Y  VECTOR  IS  (listed)" 

"ERROR  -  THE  X  VECTOR  IS  NOT  ARRANGED  PROPERLY. 

ERROR  DETECTED  BY  BICUB1 . 

THE  X  VECTOR  IS  (listed)." 

"ERROR  -  THE  PARAMETER  MAX  OF  SUB.  BICUB1  WAS  SET 
TO  _.  IT  SHOULD  BE 

"ERROR  -  THE  X  VECTOR  HAS  _POINTS  AND  THE  MINIMUM 
ALLOWED  IS  4." 

"ERROR  -  THE  Y  VECTOR  HAS  _ _ POINTS  AND  THE  MINIMUM 

ALLOWED  IS  4." 

Subroutine  BICUB2  will  print  out  one  or  more  of  the 
following  messages  if  an  attempt  is  made  to  interpolate 
a  point  beyond  the  boundaries  of  the  rtesh: 

"ERROR  -  XPT  OUT  OF  BOUNDS 

DETECTED  BY  SUB.  BICUB2" 

"ERROR  -  YPT  OUT  OF  BOUNDS 
DETECTED  BY  BICUB2" 
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3.6  Informative  Messages  to  the  User 
None . 

3.7  Input 

No  data  are  input.  See  Section  3.2. 

3 .8  Output 

(1)  Completion  of  the  P,  Q  and  S  arrays. 

(2)  The  value  of  the  function  u(x,y)  for  any  given 
(x,y)  within  the  domain. 

3.9  Formats 

Not  applicable. 

3.10  External  Routines  and  Symbols 

BICUB1  -  EDGES 
GETBP 
SOLVIT 

EDGES  *  LAGRAN 


(deck) 

ft 

II 

ft 


3.11  Timing 

The  time  required  by  subroutine  BICUB1  is  dependent  upon 
the  mesh  size  and  whether  or  not  the  normal  derivatives 
along  the  boundary  are  known.  In  the  example  of  Section 
7.0  BICUB1  took  approximately  23  milliseconds  for  a 
5x6  mesh  when  the  boundary  derivatives  were  known  and 
approximately  135  milliseconds  when  they  were  unknown. 

The  time  required  for  a  call  to  BICUB2  is  dependent  on 
the  mesh  size.  In  the  example  of  Section  7.0  an  average 
call  took  about  3  milliseconds. 


These  time  estimates  should  be  considered  very  rough 
because  of  the  method  used  to  obtain  them  and  the 
inaccuracies  of  the  timing  function  (TIMELEFT)  used. 


L SL 
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3.12  Accuracy 


An  excellent  discussion  of  the  errors  involved  in 
bicubic  spline  interpolation  is  given  in  the  paper 
by  G.  Birkhoff  and  C.  de  Boor  [8J.  Here,  we  will 
simply  mention  that  for  the  special  case  of  a 
4x4  mesh,  the  interpolated  values  u(x,y)  will  agree 
exactly,  as  they  must,  with  those  obtained  from  a 
third  order  two  dimensional  Lagrange  interpolating 
polynomial.  In  this  case,  the  remainder  term  is  well 
known  [3]  .  The  final  accuracy  is  dependent  upon  both 
discretion  and  rounding  errors.  A  rough  order  of 
magnitude  for  these  errors  may  be  obtained  from 
Section  7.0.  The  reader  is  referred  to  [4  and  5)  for 
further  discussion. 

3.13  Cautions  to  Users 


If  the  values  are  highly  variable  along  i  or  j , 

the  interpolation  "surface  may  be  forced  to  have  un¬ 
usually  high  mountains  and  deep  valleys  in  order  to 
maintain  two  continuous  derivatives  and  still  pass 
through  all  the  data  points.  In  fact,  some  inter¬ 
polated  values  of  u  may  be  so  large  or  so  small  as  to 
be  physically  unrealistic.  Whether  or  not  this  is 
the  case  will  depend  on  the  particular  problem.  The 
author  has  found  that  plotting  several  interpolated 
values,  between  and  together  with  the  given  values, 
along  a  fixed  direction  in  the  x-y  plane  is  helnful 
in  detecting  such  conditions.  In  any  event,  care 
should  be  taken  as  the  interpolation  process  could 
cause  a  physical  model  to  generate  faulty  predictions. 
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3.14  Program  Deck  Structure 


7 

g  JOB  card 
l  FTN  card 

Users  program  (containing  calls  to  BICUB1  and  BICUB2) 

Subroutine  BICUB1  ) 

Subroutine  BICUB2  / 

Subroutine  EDGES  |  E,  -  NRL  -  BICUBIC 

Subroutine  GETBP  l 

Subroutine  LAGRAN  I 

Subroutine  SOLVIT  ' 

SCOPE  card 

g  LOAD  card 

g  RUN  card 
E0F 
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4.0  METHOD  OR  ALGORITHM 
4.1  Subroutine  BICUB1 

We  begin  by  considering  the  problem  of  estimating 
the  required  boundary  derivatives  (see  Section  3.2) 
under  the  assumption  that  they  are  unknown. 

Consider  a  3rd  order  two  dimensional  Lagrange 
interpolating  polynomial  over  a  moving  (a  and  c 
are  variable)  4x4  submesh,  i.e., 


d 


Y.(y)  =  n 
J  k=c 


yj-yk 


Differentiating  Eq.  (3)  we  can  obtain  closed  form 
expressions  for  vx(x,y),  vy(x,y),  and  vxy(x,v). 

Subroutine  LAGRAN  can  evaluate  these  expressions  at 
any  mesh  point. 

Basically,  subroutine  EDGES  moves  the  4x4  submesh 
of  Eq.  (3)  along  the  boundaries  of  Figure  1  while 
calling  subroutine  LAGRAN  to  obtain  the  required 
normal  derivatives. 

Once  the  required  boundary  derivatives  are  obtained, 
the  rest  of  the  derivatives,  p.^,  q-j ,  and  S..,  are 

obtained  for  each  ij  mesh  point  by  using  the 
algorithm  described  in  [l]  ,  pages  217-218. 


14 
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4.2  Subroutine  BICUB2 


To  interpolate  a  value,  u(x,y),  at  the  point  (x,y), 
subroutine  BICUB2  begins  by  performing  a  binary 
search  to  determine  the  ij  rectangle  in  which  the 
point  lies.  A  binary  search  [7]  is  used  on  the 
assumption  that  for  most  problems  the  X  and  Y 
vectors  will  be  large  enough  to  exceed  the  break 
even  point  between  sequential  and  binary  searches 
(about  50  points) . 

Once  i  and  j  are  determined,  a  cubic  Hermite  basis 
is  used  to  evaluate  u(x,y).  That  is, 


4 

u(x,y)  *£  dr(x,h)  *s(y,k)  (4) 

r  ,s*l 

where 


“i-lj-l 

qi-l  J*1 

“i.j 

qi,M 

Si-lJ-l 

PM 

Si  J-l 

(5) 


h  •  xt-  X..J 


x  -  x 


i-1 


7m 

yM 


(6) 

(7) 

(8) 
(9) 


S 

t 
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L 


<|>3(x,h) 


4>4(x,h) 


(10) 

(ID 

(12) 

(13) 


The  functions  \p  ,  for  s  *  1  to  4 ,  are  obtained 

by  replacing  0,  x',  and  h  by  0,y',  and  k 
respectively  in  Eqs.  (10)  through  (13). 

This  procedure  requires  the  storage  of  four  values 
(u. .  p.  q..,  and  S..)  for  each  mesh  point. 

And  the  evaluation  of  u(x,y)  by  BICUB2  requires 
32  additions/subtractions  and  27  multiplications/ 
divisions.  Assuming  a  multiplication/division  is 
equivalent  in  time  to  three  additions/subtractions, 
this  results  in  113  "operations'*. 

An  alternative  approach,  not  taken  in  this  report, 
would  be  to  convert  to  a  local  power  basis. 

In  particular,  calculate  the  16  values  of 

(see  Eq.(2))  for  each  ij  rectangle,  as  described 
in  [1]  ,  and  store  them  for  each  ij  rectangle  of 
the  mesh.  This  would  require  approximately  four 
times  as  much  storage  as  the  above  method.  Also 
about  52  additions/subtractions  and  76  multiplications/ 
divisions  (270  "operations")  would  be  required  to 

obtain  the  16  coefficients  y^,  for  m,n  =  0  through  3. 
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The  advantage  of  this  approach  is  that  only  about 
19  additions/subtractions  and  15  multiplications/ 
divisions  (64  "operations")  would  be  required  by 
BICUB2  for  the  evaluation.  This  suggests  that, 
for  very  fine  mesh  evaluations,  in  which  every 
bicubic  polynomial  is  evaluated  on  the  average 
six  or  more  times,  it  is  more  efficient  to  obtain 

the  local  -power  basis  coefficients,  for  the 
r  mn 

entire  mesh  once  and  save  them.  Of  course,  this 

results  in  a  rather  severe  penalty  in  storage. 

In  contrast,  by  using  the  Hermite  basis,  as  we've 
done  here,  evaluation  costs  slightly  more  work 
but  considerably  less  storage. 

In  summary,  it  seems  best  for  most  applications 
to  save  only  the  partials  at  the  mesh  points  and 
use  the  Hermite  basis  approach. 
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5.0  SOURCE  LANGUAGE  LISTING 


(DIN?  NUMIIR  llOOilOO 
flTL|  ■  BICUBIC  SPLINE  (NTIRMlaTIIN 
(DIN?  NANI  •  U«NRL>l!CUIlC 
LANOOAQI ••PERTRAN 

cimqti*  .  eocoibo 


CBNTdtBUTBB 


JRHN  IJ  CIRNYN 

BOB  1411 

J  »•«"' 


lNf|RP«T|RN  MICIIIINQ  lYITlHl  IRANCH 
ORMMUUCAUINI  Be f Inch  oi vibibn 
IMANIIATIIN  aNRL  •  NAViL  RlHARCM  LAURA?!!? 

NAlNlNOflN,  D|Ct  tOITI 

DAT!  •  PlIRUARY  ms 


RURRtH 

IUBRbUTINII  tlouii  AND  BICUBt  RIRftRM  BtCUBtC  MUNI 
jntirpilatirn  ip  a  v aiuLatid  functibn  ir  thr  variaiu 


FUNCriBN  If  THR  VABIABLIBi 


RCPIRINcI  -  Ofi  RRIRie,  •  itcunc  MUNI  INTERRELATION 

IJ,  IP  f AYHIHATICI  AND  PNYOICIi 
VBL  41«fP  «S|*2|Ri(iR*R) 


IUIRIUTINI  CALLING  IRCIR  ••• 


THE  UIIMI  PMORAM  CALLS  lUlRlllTlNIO  BICUBt  AND/OR  |ICUItu 
SUtRlUT|N|  BICUBt  CALLS  SUBRSUTjNfi  IDOlliOETRP  AND  IRLVjt, 
IUIRRUTTNI  EDQIt  f ALLS  EUORlUTlNE  LAORAN, 

SVIRfUTfNKI  lieulS.OITIR.IlLVlT,  AnD  LAORAn  Dr  nRT  CALL 
AN?  RTHIR  JUIRSUTINII, 


SICTIBNARv  P|R  BICUBIC  SPLINE  INTIRPfLATIRN 


OICTtBNARV  CBDE  •  THE  BBCBnD  LITTiR  IP  A  LINE  |P  THE  DICTIBnARV 
indicates  the  type  rp  entry  being  DEBCRIBED. 
this  CEDE  IS  At  PILLEWI 
DICTIINARV  CSCE  TYPE  IP  VARIAILI 
CA  AERfiVUTIEN 

Cl  INTIOIR  VARIAILI 

Cl  A  INTIOIR  VARIAILI  ARRAY 

CR  REAL  VARIAILI 

CRA  RIAL  VARIAILI  ARRAY 

Cl  RRUTlNI 

CT  HERD  PRIM  TEXT 


array 


PILLRHINO  THE  DICTIINARY  CEDI  IS  A  NUMRIR  WHICH  INDICATES  THE 


C  lUIMuTiN!  NUMBER(S)  IN  WHICH  THB  ENTRY  ARREARS 

C  IU>MUT|NE  NAME  St 0 * ® U T I SE  NS , 

c  bjcuii  1 

e  B  x  cuss  2 

0  BOOES  s 

C  OITiR  4 

C  IAORAN  6 

e  SSL V I T  7 

C  CHECK  8 

C 

C  FSR  EXAMPLE,  CR2|3  KEANS  THE  ENT*v  ARREARS  IN  SUBROUTINES 

c  Btcus2  and  idoes  anc  is  a  rial  variable, 
c 

c  .  . 

970 

990 

990 

*60 

*10 

620 

630 

640 

690 

660 

670 

6S0 

690 

700 

91  A 

CSX, 6 

~ B I  CUB 1 

• 

A  SUBBRUT 1 NE  F ® R  CALCULATING  THE  N6RMAL 

7*  v 

720 

c 

DERIVATIVE!  AT  EACH  HESH  RIJNT, 

730 

CBt.S 

BICUBS 

• 

A  SUB R8UT I HE  FSB  INTERRELATING  A  VALUE 

740 

c 

LRT  AT  ANY  RIJNT  IXPT, YPT)  WITHIN  SR  IN 

790 

c 

tHI  HESH, 

760 

CRAl.4’,7 

B|Mi<I) 

9 

WELDS  THI  ELEMENT  B( 1 i I >1)  sr  EG, 19  IF  RIF  1 

770 

CRAJi4,7 

BIRK!) 

9 

HlLDS  THE  ELEMENT  BU»!*1>  |F  10,19  SF  RIF  1 

7*0 

CRAl.4‘,7 

BR<  I  > 

9 

HSLDS  THE  ELEMENT  B(!,D  |F  Ifl,  19  |F  REF  1, 

790 

C 

The  arrays  bimi.bipi,  and  sr  are  used  ti  held 

600 

C 

THE  1  AND  8  RRlME  DIFFERENCE  ARRAYS  GIVEN 

G10 

c 

In  THE  REFERENCE,  these  arrays  are  used  Fir 

020 

c 

OALSStAN  ELJMIMTISN, 

830 

css 

CHICK 

• 

A  SAKRLE  RR80RAH  TS  CHECK  BUT  THE 

840 

c 

sieuEie  SPLINE  PACKAGE 

890 

CRJ 

CXI 

9 

x  eerendent  tekrsrary  variable 

860 

CBS 

CX2 

9 

AN  X  DEPENDENT  TEMPSRARY  VARIABLE 

870 

CBS 

CX3 

9 

AN  X  DEPENDENT  TEMPBRARY  VARIABLE 

8R0 

CRS 

CYl 

9 

A  Y  DEPENDENT  TEMPBRARY  VARIaBLI 

890 

CRJ 

CY2 

• 

A  Y  DEPENDENT  TEMPBRARY  VARlABLl 

900 

CRS 

CY3 

9 

A  Y  DEPENDENT  TEMPBRARY  VARlABLl 

910 

CRA1.7 

DID 

9 

CIRRfSPINDS  T|  TUB  D  MATRIX  SF  EG,  13  SF  THE 

920 

C 

REFERENCE 

930 

CR« 

DENSh 

9 

A  'DEMHJKATER 

940 

CM 

D|Rr 

9 

DIFFENCE  BETWEEN  InTERPBLATED  AND  EXACT  VALUES 

9?0 

CRA1.7 

osm 

9 

CIRRESFINDS  TB  THB  D  PRIME  VgCTBR  BF  BO ,  16 

960 

C 

EF  THE  RIF. 

970 

CR8 

DRT 

9 

THE  ANSWER  GIVEN  BY  LAORAn 

910 

CRA1 

Dim 

9 

AN  ARRAY  FCR  HlLDlNQ  DIFFERENCES  |N  SlLVtNO 

990 

C 

BO  12  AND  14  |F  RBF , (D6BRBR) 

1000 

CR4 

DXL 

9 

delta  X  left 

ioio 

CR4 

DXR 

9 

delta  X  right 

1020 

CS1.3 

EDGES 

A  SUBREUTINE  F«B  ESTIMATING  ThI  REQUIRED 

1030 

C 

NERMAL  DERIVATIVE!  ALlNO  THE  MESH  BIUNDARIBS 

1040 

c 

ASSUHIkG  they  havb  nbt  BIEn  OIVEn.  USInG 

1090 

c 

A  THIRD  IROBR  TWB  DIMEN8IBNAL  LAGRANGE 

1060 

c 

interrelating  filtnbmial 

1070 

CS1.4 

GSTSR 

• 

A  EUBREUTJNE  FCR  CALCULATING  THE  TWB 

1060 

C 

DIHCMIIKAL  DIFFERENCE  ARRAYS  B  AND  B  RRlME 

1090 

c 

IF  EC  19  IF  REF.IDEIBBR) 

llfiO 

cis-® 

1 

9 

INDEX.  GENERALLY  USED  F|R  THE  X  ARRAY 

1110 

Cll.4,2,7 

IH1 

9 

1  MINUS  INE 

1120 

C 1 4 

CJ  jO»6»7,l 
CUi2 

eu 

cilia 


ea.i 

eu 

CS3.6 

c 

c 

c 

c 

c 

c 

c 

C l 1 *3 .  6 | 8 

c 

eiiis 

c 

cn 

CIS, 6 
C 

cu.e 

cu 

cu 

cu 

c 

CIS. 6 

c 

CIt-3,  6,8 

c 

CI7.4 

CI1.8 

C 

C 

c 

c 

c 

c 

Cl*. 8 

C 

C 

C 

C 

CIS.  A 
C 

CU.4.7,8 

Cll.4 

CU 

cii 

cu 

c 

CIS. 6 


182  •  I  »US  2 

J  •  INDEX  0ENE4ALL7  USED  F8R  THE  Y  A  MAT 

JMi  •  g  M  JNUS  INI 

JPl  •  ij  PttS 

K  >  AN  INDEX 

KH  •  I'PPfiP  LUJT  P|P  BINARY  SEARCH 

Kb  ■  LINER  LUJT  P|R  BInARY  SEARCH 

L  ■  A  CEINTER 

Li  •  L  ILLS  IKE 

LAQRAN  ■  A  SUBROUTINE  POR  DETERMINING  THE  VALUE  • t  A 
111  CIHEMUNAL  laoranoe  INTERPOLATING 
PILYMHUL  ip  arbitrary  degree  in  V  AND  Y, 

ANC  ITS  CSRI VATlVfiS  WITH  RESPECT  T«  X. 
hITH  RESPECT  T8  Y» AND  WITH  RESPECT  Ti  R|TH 
X  ANC  Y  AT  ANY  INTERSECT ! EM  P|INT  Bp  A 
V  ANDY  AT  ANY  INTERSECTUN  PBINT  BP  A 
THI  CUEMiBNAL  MESH, 

M  •  THE  NUMBER  IP  Y  MINTS  AT  UHtCH  THE  PUNCTIBN 

NAS  8BSERVED.  MysT  BE  GREATER  THAN  3, 

MAX  •  THI  CREATE"  IP  N  AND  M, 

HAXS  ■  THE  VALUE  MAX  SHflULD  HAVE  BEEN  SET  T« . 

HP  •  NUMBER  IP  THE  PINAL  MINT  BN  Y  AXIS  TI 

BE  USEC  }N  LAORANQE  INTERPBLATIBN 
MM1  ■  M  MINUS  BNE 

MM2  •  M  M  JNUS  TWB 

MMlN  •  THE  SMALLEST  VALUE  M  CAN  TAKE, 

MPT  •  THE  VALUE  BP  THE  Y  VECTBR  Tfl  BE  USED  IN  LAORANOE 

INTERPILATIIN 

MS  •  NUMBER  BP  THE  STARTINQ  PBINT  BN  V  AXIS  TO 
BE  UREC  IN  LAQRANOE  INTERPOLATION 
N  •  THE  NUMBER  |P  X  PIInTI  AT  WHICH  THE  PUNiCtlON 

NAS  EBSERVED, 

N  ■  NUMBER  bp  elements  in  LINEAR  SYSTEM 

NDPAULT  •  A  PARAMETER  WHICH  MUST  BE  SET  TB  1  IP  SUBROUTINE 
8ICUB1  1!  T|  c*Ll  EDGES  T|  CALCULATE  THE 
RECUIRED  NORMAL  DERIVATIVES  ALBNO  THE  BOUNDARIES 
BP  THE  MESH,  ]P  NDPAULT  IS  NIT  SIT  TB  1,  BtCUBl 
ASSUMES  THE  NORMAL  DERIVATIVES  PON  THE 
BOUNDARIES  HAVE  ALREADY  BEEN  ENTERED  INTO 
ARRAYS  P.Oi  AND  S  BY  THE  USERS  CALLINQ  PRIOR AM , 
NfRRBR  •  AN  ERROR  INDICATOR,  IP  THE  POINT  (XPT.VPt) 

C«l$  NIT  LIE  HITHIn  THE  MESH.  nERROR  hill 
BE  SET  TI  1  BY  8ICUB2  OTHERWISE  IT  WILL 
RBMAIN  SIT  TB  0,  IP  NERROR  {S  RETURNED  AS  i. 

AN  INTERPOLATED  VALUE  IS  NOT  COMPUTED, 

NP  a  NUMBER  |P  THE  PINAL  POINT  pN  THE  X  AXIS  t« 

BE  USED  IN  LAQRANOE  INTERPOLATION, 

NM1  •  N  MINUS  ONE 

NM*  •  N  MINUS  TWl 

NMj  •  N  MINUS  S 

NM I N  •  THE  SMALLEST  VALUi  N  CAN  TAKE 

NPT  •  THE  VALUE  IP  THE  X  VECTOR  TO  Bl  USED  IN  LAORANOE 

INTERPOLATION,  XPTaX(NPT) 

NS  •  NUMBER  IP  THE  STARTING  POINT  IN  T HE  X  AX|S  TI 
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jlBO 

1140 

iito 
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1200 
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1240 
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1*80 

1270 
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1900 

1970 
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1640 
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1670 

16R0 
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eu 

c 

c 

c 

c 

c 

c 
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•E  USEE  IN  L*QR*N8E  INTERPOLATION, 

SET  TS  1  IF  LAQRAN  IS  T0  INTERP®LATB  A 
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SET  Tl  2 

set  Tl  3 

SET  Tl  4 


CRAlO’i 
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P 

C 

CR2 

PJ1J1 

CR2 

FIJI 

CRAl-Ji 

Oil 

0 
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R 

e 

CS8 

RANF 

c 

c 

CRI  , 

RJNV 

CRAlO< 

M 

S 

C 

CStil 

SOLVIT 

C 

C 

CRA1 

ITIMP 

CRB 

BUM 

CR4 

SUMX 

CRB 

SUMY 

CRB 

T 

CR2 

71 

CRB 

TJMt 

CSS 

TJMELEI 

C 

c 

CRI 

TRJ 

CRI 

tpimi 

CRAlOi 

Oil 

U 

C 

CRI 

uixact 

C 
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UI1J1 

CR| 

UIJ1 

CRI 

UjNT 

C 

CRI 

URT 

CRA4 

V 

CRI 

Ml 

CRI 

Hi 

CR0 

HJ 

CRI 

• 
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• 
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H? 

CRA1-9, 

Bil 

X 

T(  QET  PARTIAL  DIRlVATtVE  U/R  tl  tt 
T|  G6T  PARTIAL  DERIVATIVE  H/R  Tl  V, 
T|  GET  PARTIAL  DIRIVATIVI  H/P  Tl 
lltH  X  AND  V , 

THI  NSRPAL  DERtVATIVPS  |F  U  WITH  RESPECT  t|  X 


P(I»1»J*1) 

P(!,ij«l) 

THE  NIRPAL  DERIVATIVES  IF  U  W| tH  RESPECT  Tl  Y, 
A  TEMPORARY  VARIABLE  USED  In  FIRMInQ 
TM|  C  PATRIX 

A  COC  3100  RANDOM  NUMBER  QINERATBR 
GENERATES  UNIFORMLY  DISTRIBUTED  RANDSM 
NUMBERS  BETWEEN  0  ANDi 
THE  INVERSE  OF  R 

TM|  NtRMAL  DERIVATIVES  IF  U  WITH  RESPECT  T0 
MTH  X  AND  Y, 

A  SUBROUTINE  FSR  SOLVING  A  LINEAR  SYSTEM 
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MAXIMUM  LENGTH  IS  DETERMINED  BY  THE  AMOUNT 


16«0 

i?io 

1710 

17*0 

1730 

{74o 

l7So 

£740 

17t0 

17R0 

17*0 

iooo 

1810 

{•90 

18S0 

1040 

1890 

{•BO 

1070 

1000 

{890 

1900 

1910 

1910 

{930 

£940 

1990 

1940 

1970 

1480 

1990 

9000 

9010 

9010 

9030 

9040 

9090 

9040 

9070 


9080 

9090 

9100 

9110 

9190 

9130 

9140 

9190 

9190 

0170 

9180 

9190 

9200 

9210 

9290 

9230 

2240 


21 


c 

EF  C6RI  AVAILABLE! 

2250 

CRi 

XJI 

• 

X(3)KJSUS  »(2) 

9260 

CR6«6 

X| 

t 

THE  ITH  VALUE  IF  X 

92*0 

CR| 

X I  Hi 

a 

X(M> 

22*0 

CRI 

XMXlflH 

• 

TlHRgRARV  VAR1A8LS 

22*0 

CRI 

XN12 

• 

X  (NH1 >«X(NH2 ) 

>300 

CRB 

X> 

* 

XURT) 

9310 

CR2 

C 

XRT 

• 

THE  X  CI»I*D  INATI  §F  THE  RBINT  T8  BE 
INTERRILATED, 

>320 

9390 

CRAJ«3. 

C 

C 

C 

C 

6*1  T 

• 

THE  VSCTCR  IF  DISTINCT  VALUES  IF  THE  SECOND 
INCEPENDf NT  VARIABLE  ARRANOED  IN  ASCENDING 
CRCER.THE  minimum  length  IF  Y  IS  4 

AND  THE  MAXIMUM  LENQTH  IS  DETERMINED  BY 

THE  AMOUNT  IF  C|Rl  AVAILABLE. 

9340 

2350 

>360 

>370 

>360 

CRi 

Til 

V 

YCJMI2) 

>3*0 

CRB 

Tl 

9 

THE  ITH  VALUE  SF  Y 

>460 

CR6 

TJ 

■ 

y  (y) 

>410 

CRI 

TJMl 

a 

Y  Cj-l ) 

2420 

CRI 

YM 12 

a 

Y<HMl>-YiMH2> 

>490 

CRI 

VMYlBK 

a 

TlMRgRARV  VARI AILS 

1440 

CRB 

T> 

• 

VIMRT) 

>450 

CR| 

C 

V>T 

a  THf  V  ClafROlMTE  BF  THE  MINT  TB  BE 

intemiiatid, 

>460 

>470 

CRA7 

e 

Z 

a  THE  S6LUTI6N  VECTOR  FIR  THE  LINEAR  SYSTEM 

>410 

>4*0 

>560 

•UMIUTfNK  BtCUBKK  i>>  i  N  Of  AULT ,  X  ( V » U  i  ^ « 0 1 1 « HAK .  BP « B  t  PI  •  •  1  HI*.  O’.  Dl  • 
X  IT|m>,dR) 

BICUBIC  IBlINl  InTERRILATJIN 

THIS  SUBROUTINE  C A L CUC A T( 3  THE  MARTIAL  DERIVATIVES  FIR  THE  HISH 

DIH|N|IbN  X(K)iV(M)«X(K»T<)*tPCN«H>iO(M|N};S(HiN) 

D I M|N$ 1 0H  BR(HAX).B1*1(MAX).B1M1<MAX)(D(MAX),DS<MAX),STEm><MAX)i 
X  DHfMAX) 

0ATa<NM|N|4),(HMIN«4) 

CHECK  T|  SEE  |I  THE  max  EARAKETER  WAS  SfiT  CIRRICTlY 


MAXS  IS  HHAf  MAX  SHBULt  IE 
MAXfaN 

IF(N,IT;m)  MAXSaM 
IFIMAXS  , NS ,  HAX)  66  Tt  «00 


DETERMINE  WHETHER  THE  X  ANE  V  VECT8RS  aHE  WITHIN  LIMITS 


ir (N  >LTt  WHIN)  OB  Tf  BOS 
IF<H  |LT|  MHIN)  06  T|  BQT 
NMliN.i 


>510 

isto 

ESSO 

>540 

>550 

>560 

>5T0 

>560 

>5*0 

9650 

?610 

9620 

2690 

>640 

9650 

9660 

96T0 

>680 

>6*0 

>700 

2710 

9720 

9790 

9740 

9750 

9760 

2770 

>7*0 

27*0 

9600 
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o  <r»c»  o»  non  nnnn  n  non  non  nnnn  non  n  nnnn 


NH2|Na2 

HHlfHal 

NNllHit 

DETERH tNE  WHETHER  THE  X  INC  Y  ARRAYS  C8NYA I N  DISTINCT  ELEMENTS  AND 
ARE  ARRANQlD  IN  AtCINOINQ  IRCER, 

Dl  100  |«1»NM1 

IF (XCI )  .OE,  X( 1*1) )  01  TS  911 
too  CENT  1 NUE 

01  tlO  JiliHHl 

ir<V(J)  ,0E ,  Y(J*D)  01  TS  912 
tlO  C9NTINUE 

determine  the  eooe  boundaries  fir  Miand  i  i f  requested 
it (ndfault  ,eo,  d  call  edoes(n,m,x,y,u,r,o,s) 

OET  THE  DIFFERENCE  aRRa Y>S«  FIR  THE  X  VECTOR  AND  ALSS 
OET  THE  B  PRIME  ARRaY.IR,  FRIM  TNI  I  ARRAY 

CALL  OSTBFIN.X.lR.lIRldlHi) 

NSW  SOLVE  F|R  THE  PAROLS  W/R  T<  X  WHICH  WERE  NST  QtVlN 

X32|X|3)*X(2) 

XM2*X(NHD«X(NH2) 

Dl  30  Jal«M 

SET  UR  THE  D  VECTOR  F|R  10(11)  |F  R|F 

Dl  SO  I  Hill  a  NH2 

IllMlAl 

IF1|I*1 

R«(X<|)|X(tHl))/(X(IRl)<X(i)) 

RJNV*1,/R 

01  iHl )*3i*IR*<UC  IRiiiJ>«UC  1 1 J)  )*RlNV*(UI  (’•  J)«U(  (Mil  J) ) ) 

IF(J  ,10.  1  ,IR,  J,BC.  M)  34,g9 

34  DS(!Ht)93,*(R*l0lj',IFt>«C('j,m*RlNV*(0lJ,n.0IJ#INin) 

note  0  AND  S  ARRAYS  ARE  (TIRED  AS  0(U,I)  AND  S(J,I)  RATHER  THAN  , 
0(I|J)  SECAOSI  IF  FORTRAN  CONVENTIONS  FIR  (TORINO  ARRAYS  COLUMNWISE 

39  CINTINUl 

ADD  ADDITIONAL  TERMS  T|  THE  FIRST  AND  LAST  ELEMENTS  IF  THE  D  VSCtlR 

D(1)*D<1>»X32*R( tAj) 

D(Nh2)*D(NM2)«XN12*R(N,J) 


NIW  SSLVE  LINEAR  SYSTEMS  FIR  E0(H)  OF  REFERENCE, 
CALL  SH.VlT(NN2,R<2,g>,D,M,flIFliBJMiiDF) 


mo 

2020 

2630 

2140 

2190 

2140 

2020 

28*0 

21*0 

2900 

2910 

2920 

2930 

2940 

2990 

29*0 

2920 

2960 

2990 

3000 

3010 

3020 

3030 

3040 

3090 

30«0 

3020 

30*0 

30*0 

3100 

3110 

3120 

3130 

3140 

3190 

3140 

Sl20 


3210 

3220 

3230 

3240 

3290 

3240 

3220 

3260 

3290 

3300 

3310 

3320 

3330 

33*0 

3390 

3340 


r 


ir C J  ,l0,  1  ,11),  J  Jo,  M)  37,30 

e  too  additional  tiohi  re  the  os  array 

c 

37  Dt<l>a0s<l).X32*Stj,}) 

OSCNH3  )«Of <NH2).XNi2»9 ( J,N  > 

C 

C  NIW  SOLVE  LINEAR  SVITBP  MR  EQ(l2) 

C 

CALL  MLV|T(NH2ilT2HR,CS,|R,B!Rl,BfMl,DR) 

C 

C  HIVE  VALUE!  MBm  The  TIhMRARV  ARRAY  InTS  The  RARTiAL 
C  ARRAY  f|R  The  CR8SS  TERMS,  HE  Cl  THIS  BECAUSE  07  F0RTRAN  ARRAY 
C  ST0RAO6  CONyONT IONS 
C 

DE  10  Lal«NN2 
LHL*l 

l<J,Li>*»TlHMLI 
10  CINTINUI 
C 

30  CBNTINUE 

C  .  , 

C  BET  THE  DIFFERENCE  ARRAY  ,1,  FAR  THE  Y  VECTBR  AND  THE  I  PRIME  ARRAY 
C 

CALL  QI?BP<M,V,BP,B!P1,8!M1> 

C 

C  NEW  OET  THE  PARTIAL!  HiTH  RESPECT  TS  Y 

e  WHICH  H|R|  NET  GIVEN,  I,  E,  MEMBERS  SF  THE  G  ARRAY  HHIOH  HIRE  Nit 

C  specified 

C 

Y32«V(3)*V(2) 

YMia«Y(MMl)»Y(MM2) 

DE  40  Ial|N 
C 

C  SET  UP  THE  D  MATRIX  F|R  |Q(13)  |F  REF 
C 

DE  41  JH1B1,MM2 

•jajMiai 

jPiaJ«l 

Ra(Y(  JlfY(  JM1M/(Y(  JPl)«YtJ)  > 

«INV»1./R  , 

0UMiM3,a<R*<U(!iuPi>>UUH>>«RlNV«(U(!,J>*U(!'JMim 
DS(jMi)a3,«(R*<Pii,jPi)*p(t,gn»RiNVa(p(i,j)>p(r,jMim 
49  CENTINUI 
C 

D(l)aO(i>aVJ2*S(l«n 

D(MMl)aDJMME».VMl|aO(M,n 

C 

0t(l)SDf(l)-Y32*S(|,n 

D8(HH2>*Dt<HM2).YHl2*S(M, I ) 

C 

C  NEH  CAN  SILVB  THE  LINEAR  EyITEM  EF  «0(lJ>  PER  THIS  i 

c 

CALL  fELViT(MM2,0(2il),D,IP,EIPlilIMl,DP) 

C 


3370 

33B0 

33(0 

3400 

3410 

3410 

3430 

3440 

3490 

3440 

3470 

34*0 

3400 

3900 

3910 

3920 

3930 

3940 

3990 

3940 

3970 

39R0 

3900 

3490 

3410 

3420 

3430 

3440 

3690 

3640 

3470 

3600 

3600 

3790 

3710 

3710 

3730 

3740 

3790 

3760 

3770 

37Bb 

3700 

3E0O 

3010 

3SIQ 

3E30 

3E4Q 

3E90 

3160 

SOTO 

3110 

3B(Q 

3090 

3010 

3020 
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o  ao 


e  \sw  stive  the  linear  systekssf  eg<i4>  fir  this  i 

c 

CALL  SSlYJTlMMJ.Slj.n.DS.BR.BlRl.SlMi.DR) 

40  CSNfINUE 
C 

c 

RETURN 

c 

100  RRInT  9g3|HAXiMAXS 
STAR 

•OS  RRINT  9o6|N,NH|N 

star 

•07  RRInT  9oMiHh!n 
'  STAR 

•n  rrint  •is,mn,l»itM 

STAR 

•12  RRINT  «l6|(Y(U)iJ«ltW) 

STAR 

«•••* «”t»t*iF|RMAT  STATEMENTS  ••••*«••»•••»**•••••••••••»••**•• 

•03  FSRMAT(10X,«ERR|R.TH|  PARAMETER  max  BF  Sub,  BICUB1  HAS  SIT  Tl  •, 
X  1 9 , * )  [T  SHAULO  BE  *i IS) 

•06  FSRmAT(10X.«|RR|R-THI  X  VECTOR  HAS  MSl 
X  *R|lNTS  AND  THE  MINIMUM  ALLSHBD  IS  *iI3) 

90S  F|RMAT(10X(»§RR|R*TH|  y  VECTSR  has  MS. 

X  *R|INTs  AND  The  MINIMUM  ALLSHBD  II  •* IS) 

•IS  rsAMAT<l0X.«SRR0R«THf  X  VECTSR  IS  N*T  ARRAnOID  RR0RERLV*',//, 

X  iox,«irrsr  DETECTED  BY  BICUll*./, 

X  IOX.RThE  X  VECTSR  IS*./. 

X(S<2X|Fi2.4j  J ) 

•16  F|RMAT(10X.*IRRSR.THE  Y  VECTSR  H  NIT  ARRANOID  RR0RERIY*.//, 

X  10x.*ErR8R  DETECTED  BY  SlCUBl*,/, 

X  10X.*Th8  Y  VECT9R  IS*./. 

X(S<2X|Fi2.A) ) ) 

C 

BNO 
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aon  ooo  c*  a  anna  nn 


fUBMUTjNK  l!CUl2tXAt,VAY,'l'Af  »N|AA|A',NiM»X»Y,U»A»0»I) 

•UlMUTtNK  yl  iNTlApHATI  I  VAIUMaT,  AT  LlCATJlN  <«pT.YpT» , 

THIS  IUIMUTINB  AligMtf  THAT  THl  X  A*0  Y  VICTIM  aA|  LAAOI  INIUQM 
tl  HAKI  A  lfNA*r  IIA»0K  YICHN I OUl  lU’IAUA  Tl  A 
9I0UBN T 1 AU  8iA»CH  YICHNtOUf, 

N|AA|A  j»  IfT  T#  1  jr  (XATjVAT)  tl  lUTllDI  Th|  AA|A|A  D|H  A I N 
D ImInII |N 
NIAMAao 


tlNDUCT  A  IlNAAy  «lA*0  t|A  t 

KH«N 

Kt«l 

2  tl(KLAKM)/* 

3  tr(xd).XAT)  4,T,li 
11  IKXAT-XM-U)  4.1M 

4  KHA| 

9  rr (khakl'-d  ««t<> 

6  KLA| 

01  Tl  9 

9  tr<XAT-X<KM>>  10,1,13 
8  tAKH 
01  Tl  T 

10  tr(XAT-X«KD)  1 3 1 14 1 1 
14  IAKL 
01  tl  7 
18  IaIaI 
01  Tl  7 
13  NIAAIAai 
PAINT  12 

(JAt|«0.0 

AftgAN 

eiNDUCT  a  iinaay  siaaoh  ma  g 

7  KHAH 
KL«t 

20  J«<KlAK*)/2 
30  |P<TCJ)iVPT>  IO.YO.liO 
iio  fr «yPT»Y« J-ll )  40.110, TO 
40  KHAJ 

90  tnxHAKL«|)l6, 90,20 
10  KbAJ 

01  tl  lo 

90  tT(YAT«T(KH))  100,10,130 

"  «“?.  T, 

100  ir<YPT-Y<KL» »  130,140,10 
140  JAK|, 

01  Tl  TO 
ilO  jAJtl 
01  tl  TO 
130  NIAAIAai 
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nonan  ci  nnoacioo 


mint  13 
UPTmOiO 
MTuRN 


uii  a  egiie  mhniti  iaih  ti  bvaiuati  Thk  bicuiic  piiynomjh, 

IT  (XPT.YPT)  ,  .  . 

(XM.YM)  UII  NITNIK  THI  MCTANOll  BIUNBID  IV  X( 1 1 .X< Ml ) . YI J). 
AND  Y(J.l), 

70  I H 1 • 1 ■ 1 
JHlfJ.l 
X|»XCJ) 
yj»y<j> 

X(MiaX(lH|) 

jUj) 

XHX4«H*|X|T.X|Ni)/(X}«XlHl) 

YKYi«K«( YPT.VJM! 1/ ( Yj.YJNl ) 

CXli(3,i2,*XHXl|H)iXKXllH**a 

CYl«(J,iJ,«YMYHK)»Y»»YilX*»j 

CXIt<XliXRT»MNXi|N 

eY2<( VJ|YIT)*VNVllK 

exsiXMx«iH*ex2 

CY3(YNYjOK*CY2 

eX2|CXItCX3 

CY2aCY2aCY3 

UUl«U<T»JMl) 

RUlM«r»JMU 

R|lJI«R|IN&tJN|| 

TIMiagiiJi*CYlMUUhi^».UIiJl>*CY|*0IJHl,lHl)»cy3«0(J.lHl) 

TJayl J1*CY1*(U( I t J)  -ttJl>*CY2*0(jHiiI).CY3«0(J,n 
T7lHllRllJl*CYl*(^(lKl,J>.PUJi)*CY2*l(jHi,IHl).CY3*S(J(!Hl) 
YTIlRIJt*ri»«T(IO>.MJl)lCY2*|(JMl,I)»CY3«l<J,n 
UITaYINj«CrXt4{YJ«TlMl  >*CX2»TPlMt"CXl*TPl 

MTgRN 


PMMAT  statiminti 

. . 

12  MRmAT(iOXi«MMX*XPT  CUT  I F  8IU^DI*i/i 
X  10Xi*O|T|CTI0  IV  SUI,  Itoullo 
23  fMMAT(tOXi«MR|R*YPT  IUT  IP  IIUNDMi/'i 

x  iox««oiTseTiD  iy  ireuia*) 
c 

180 


AMO 

4170 

AMO 


AMO 

AliO 

AI10 

4110 

4130 

AMO 

4110 

4100 

4170 

AMO 


4110 

»oio 

*010 

*oio 

*030 

mi 

9040 

9070 

9010 

9010 

9160 

*110 

9110 

9130 

9140 

9190 

9140 

9170 

9110 

9110 

9260 

9210 

9210 

9230 

9240 

9290 

9140 

9270 


9210 

9210 

93d0 

9310 

9310 

9330 
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L. 


o  on  aoo  onoaann  aa  nn  n  r>  a  nci 


SUBRBgTjNl  EDOEIIN',H,XiYiU>iO.I) 

tr  THE  USE*  MOUIITID  DIF4 ULT  SCO!  C«NDtTt«NS 

THIS  JUIMUTINI  WILL  DETERMINE  THEM 

US  I  NO  A  LAOrANOE  INTERPOLATING  PlLYNlMlAL  8F  ORDER  3X3 


THIS  lUlMUTlNE  IS  GALLIC  By  lUlMUTlNf  BICUR1, 

THIS  SUBROUTINE  CALLS  SUSNCUT J NE  LAQRAN, 

DlMENlIgN  X(NM<M)iL(N|MhP<NiN)|0(H|N)',S(H»N> 


SET  PARTIAL!  WITH  RESPECT  TB  X  ALBNQ  EDGES 


DETERMINE  TmE  LBCATIBN  IP  THE  4  X  4  GRID  F8R  THE  LAORANQE  POLYNOMIAL 
DS  10  J«1,M 

IFImB  ,oT|  ( M*3 ) )  M|iM*3 
IF (MS  ,lT|  d  ms«i 
MF*MSp3 

CALL  LAflRAN(3il,4|><S'iMfiR(i#Vlil»J/N|M^XjY#U»RiO;SJ 
CALL  LAoRAN(2iN«3iN#MSiMF|R(N#J)#NiJiNiM,XiY#UiRiQ»S) 

10  CINT1NUB 

BET  FART  I ALS  WITH  RESPECT  TB  Y  ALBNQ  EDOES 
Dl  10  Ill.N 

nsri*i 

IF(Nl  iQT | (N*3) )  NSRN-3 
!F<nS  .LT,  1»  N8«l 
NF«nS*3 

CALL  LAQRAN(3#NSiNP»li4iC(i  #  l)tl<l/NiM/X>YiUiP«QjiS) 

CALL  LAQRAN(3iNSiNP»H«s»PiC(MiI)»l#M»N»M|XiYiUrP»Q»S) 

20  CONTINUE 

SET  PARTIAL!  WITH  RESPECT  TB  X  AND  y  AT  CBRNIRB 

CALL  LAaRANtAa.A.iiAiSIlijhiltliN'.M.X^UiPieiS) 

CALL  LAQRAN(4|Na3iNil«4iS(l|N)«NiliW»MiX>YiU»Pi0iSI 
CALL  LAQRAN(4,ll4,Hi3<K,S<»',l>(l,M^N,h,X/y(U,P.Q,8I 
CALL  LAqRAN(4«n*3iNiM*3iP»S(MiN) iNiM*NtMiX« ViUiPiQiS) 

RETURN 

END 


340 

390 

340 

370 

3*0 

3*0 

400 

410 

4«0 

430 

440 

490 

440 

470 

4*0 

4*0 

900 

910 

930 

930 

940 

990 

940 

970 

9B0 

9*0 

480 

410 

430 

410 

440 

490 

440 

47o 

4*0 

4*0 

700 

710 

780 

730 

740 

790 

740 

770 

7*0 

7*0 
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on  n  o  o  o  oo  oo  oo  o  ooo  oo  oo 


lUMIUTTNf  QITIR<N,V‘,IM!Rl«liMl> 

wimutjnI  it  ost  thi  t  iimekiiikai  oirrERENCS  array.i,  »r  socisi  «r 

RiriRENcE  AND  AlM  THE  |  MINI  ARRAY  MR  8AUISUN  IlIHJNATJIN 
tRtl  IUIRIUTINR  !!  C*LllC  IV  IDlRfUTlNI  HCUE1, 


fl  SAVE  ITERAOE  WE 

eteri  element  § < ! , i >  in  ikd 

ITEM  ELEMENT  IK  1161(1) 

iters  element  iiiiM)  u  iimkp 

THIS  REOUCH  THE  MACE  RIQIIRED  MR  THE  R  AND  •  PRIME  ARRAYl 
RREM  2  *  N  *  *  2  T|  J*N  W»RtS 

OlMlNtllN  BM«N).BlMl(N)»IIwl<N)|V(N> 

NMl«N.l 

NMEgNtE 

NMSrNbS 

DXR|V(2)«V(1) 

•R(i>Ei;*(DxR*iv(S)tv(im 

I|»1M>|DXR 

BXL|V(NI*VINMI) 

HMKKMjliDxi 

IM(NM2>|2,»<DXL*<V(NW1>«V(KM2))) 

IP CN  ,E0.  4)  BE  T|  n 

DE  10  I«2|NMS 
lMlflal 
IRlllAl 
IM2|f*l*l 
DXL«V(  1|»2|»VC  IP1> 

DXR|V(1r1»#VU> 

ItMilDiDXL 
BP  C 1 )I2|*(DXL*DXR> 
i|Ri(|)|BXR 
10  CONTINUE 


NEW  determine  THE  I  PRIME  MATRIX 

11  01  20  I|2|NM2 
IMlilal 

BPtpflRdNIlMKlHIiPKIMlI/lpIlMl) 
20  OENflNUE 
RETURN 
END 


•do 

010 

110 

190 

140 

190 

M0 

•70 

MO 

8*0 

900 

910 

910 

990 

940 

990 

960 

970 

960 

990 

odo 

010 

010 

090 

0«0 

090 

040 

070 

040 

090 

ido 

no 

120 

190 

140 

190 

1*0 

1T0 

160 

1*0 

2(|0 

210 

220 

290 


>240 

>290 

>260 

270 
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ri  ci  r»  r»  o  non  n  nnnnnnnnnnnnnnnnnnnnnnnnnn 


SUBROUTINE  lAGRAN<nTYPEiNSiNF,MS#MF»DPT»NPT,MPT,NiM,X,Y,UiOiO«B)  *2*0 

*3*0 

*5oo 

SUBROUTINE  f 00  INTERPOLATING  THE  VALUE  IP  A  FUNCTION  ANQ  |TI  *310 

DERIVATIVE!  WITH  RISPfCT  Tf  X,  WITH  *E8PFCT  T|  V,  AND  WITH  RESPECT  Tl  6320 
X  AND  V  AT  THE  MESH  MINTS  USJNO  A  L*0RANQE  INTERPOLATING  *330 

MLYnRmJAL  ®r  ARBITRARY  f  RCEP  *340 

*310 

THIS  SUBROUTINE  IS  CAUSE  BY  SUBROUTINE  EDOES,  *3*0 

*370 

*3*0 

ns*  number  er  startjno  pijnt  along  the  x  xa j s  *300 

NF*  r INAL  print  along  the  X  AXIS  *400 

MS*  ST  ART  J  Nq  print  AlINQ  THE  Y  AXIS  *410 

mf*  final  print  along  the  y  axis  «4to 

INTERPBLATIRN  is  CARRIED  ELY  BVE*  THESE  PUNTS  *430 

*440 

*490 

subroutine  laqran  has  i6  parameters  In  its  calling  sequence  *4«o 

*470 

*4*0 

NTyPB»l  get  function  ITSELF  *490 

NTVPE*2  OET  PARTIAL  WITH  RESPECT  TO  X  *900 

NTYPE  *3  GET  PARTIAL  WITH  RESPECT  TO  Y  *910 

NTYPEI4  OET  PARTIAL  WITH  RESPECT  TO  X  AND  Y  *920 

*530 

*9*0 

DIMENSION  XlNI.YlMhUHlWI.PIN.HliQlMtNl'.SlM.N)  *990 

*9*0 

XP*X(NPT)  *970 

YP*Y(MPt)  *9*0 

DPTtO,  *9*0 

00  TRU0i20,30. 401, NTYPE  *600 

*610 

CALCULATE  The  VALUE  SF  The  FUNCTION  AT  the  mesh  print  JNPT.MPT)  *«>o 

**30 

10  DPT«U(NpT|MPT)  *640 

RETURN  **90 

6**0 

*670 

CALCULATE  ThE  FIRST  PARTIAL  WITH  RESPECT  TO  X  ««*0 

AT  THE  MESH  PRINT  (NPT.MPT)  «**0 

*700 

20  DR  211  I*NS«nF  *710 

SUMaO.  *720 

DO  219  L*NS|NF  *730 

IFtL.EO;  I)  OR  T0  219  *740 

P*RD*i«  *790 

DO  212  k*NS,NF  *7*0 

IFIK.EO;  I  ,0R,  K.EQ,  L  )  00  TO  212  *770 

PR§0*PRrDRIXP»X(H) )  *7*0 

212  CONTINUE  *7*0 

SUM|9UM*PR9D  ARdO 

219  CONTINUE  *010 

DENRMlt;  4R20 

X I «X  < 1 >  *030 
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C«  219  K'NS.NF 

0040 

IF (k  ,10,  I)  08  T|  213 

8090 

DEnjHiDEN|H*<XM(iO) 

8000 

213  CONTINUE 

80TO 

BPT|DPT*8UM/0EN0H*L(liPPT) 

8080 

2 1 1  CONTINUE 

8090 

return 

8900 

8910 

CALCULATE  The  F  JRST  PARTIAL  WITH  RilPlCT 

TO  Y 

8920 

AT  T«E  HfiJH  PUNT  IN^T.HPT) 

8930 

8940 

30  D|  311  J«hS,mF 

8990 

8UM|0, 

01  319  LbMS.mF 

8980 

89T0 

IF (l  (EO.  J)  G8  T»  319 

8980 

PROD*!. 

8990 

08  312  k*HS|MF 

7060 

IFCK.SO;  J  |8P,  K ,10,  L  )  68  TO  312 

7010 

PR8D«PPeO»(YA»Y(K)> 

7020 

312  CONTINUE 

7030 

9LM«SUH*PR9D 

7040 

319  CONTINUE 

7090 

OENOHii; 

7080 

VJ«V<J) 

7|)70 

DO  313  k»MS.*F 

70"0 

IF(K,60;  J>  00  Tfl  313 
DEN0HlDEN0H«(YJ.T(Kn 

7090 

7100 

313  CONTINUE 

7110 

DPTtnPT«$UN/DENOM«L(NPT|.J) 

7120 

311  CONTINUE 

7130 

RETURN 

7140 

7190 

CALCULATE  the  PARTIAL  Min  RESPECT  T|  X 

AND  Y 

7180 

AT  THE  ME5M  POINT  (NPT'.MPT) 

7170 

7180 

40  DO  41  IiNIiNP 

7190 

sumxrO. 

00  43  LaNtiNP 

7200 

7210 

ifil.eo;  n  no  to  43 

*220 

PRfl0«l. 

7230 

DO  44  KftNI'NF 

7240 

ifik.eq;  1  ,qr.  mb,  l>  Qo  to  44 

7290 

PR8DpPR0Dp<XP-X(K)) 

7280 

44  CONTINUE 

7270 

3UHXR3UMXPPR80 

7280 

43  CONTINUE 

OENflHli; 

*290 

7360 

XtRXd) 

7310 

DO  49  K|N|,NF 

T320 

ir <K  t lQ •  I )  00  TO  49 

7330 

OE J0MIDEN0M«(XI.X(KI) 

7340 

49  eO^INUl 

7390 

9UHX«|UMX/DENOH 

7380 

DC  42  J«H||HP 

7370 

SUM'iOi 

7380 

08  433  LRM9*hT 

7390 
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!P<liiO;  J)  01  TO  433 
• 

01  434  K«M3,Mf 

|F(K.IO;  J,**.  K.BCi  L)  01  Tl  434 
RR*DMR0D*<Y*»V(*O) 

434  C|NT!NUB 

lbMyi|UHV*M«D 
433  CONTINUE 
OENOHlt; 

VJ«Y< J> 

DO  439  K'Hl.Hf 

IF (K  ,EQi  J>  <S«  PI  439 

DEN|HlOEN|H*(YJ«V(K)) 

439  CINTINUp 

9tHy»|UHT/DEN0H 

0PT*0PT49UHy*SUMX*lU»'J> 

42  CONTINUE 
41  CONTINUE 
RETURN 

END 


400 

410 

410 

430 

440 

490 

4*0 

4?0 

4*0 

4*0 

900 

910 

920 

930 

940 

990 

9*0 

9T0 

9*0 

9*0 


*00 
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X  ^ 


subroutine  ietviT(K,!(c»iMJPi,iiHi,oP)  mo 

Tt!0 

7*10 

THIS  RigTINC  |l  CAUID  ly  SUBROUTINE  ItCUIl  7*40 

7*90 

DlH|NiJ#N  BR(N),BIBl(N)iBlNl(N)|D(N)iZ(N),DR(N)  7**0 

7*70 

fcOMPUTB  TNI  D  RR1HB  VECTOR,  IEI  EQ(1*>  BP  REP  7**0 

7*70 

DP<l>lD(l)  7700 

D9  10  !«2«N  7710 

IHlllal  7790 

|R(|  )|D| I )«I|M1(  I  )*0P (!N1)/8P( IR1)  7730 

10  CONTINUE  7740 

7790 

•BUIS  HLUTI0N  By  REClRS  JIN  RIL*T!0N  OP  E0(l7)  ,  SEE  REP  77*0 

7770 

!(N)aDR(N)/BR(N)  77R0 

NHliNal  7700 

DO  20  J«1,NM1  7*00 

taNiJ  7010 

IPllUl  7020 

Z( I )*(Dp( |  )«BIRi<  1  )*Z( IRl) )/IP( t )  7830 

20  CONTINUE  7840 

RETURN  7890 

END  7840 


6.0  COMPARISONS 


For  most  quantities  defined  over  a  two  dimensional 
mesh,  where  N  and  M  are  greater  than  4,  the  bicubic 
spline  generates  a  more  physically  plausible  inter¬ 
polation  surface  than  a  two  dimensional  Lagrange 
interpolation  polynomial  over  the  same  mesh.  When 
N  and  M  areeaual  to  4 the  bicubic  spline  and  Lagrange 
interpolating  polynomial  are  identical . 

No  comparisons  have  been  made  with  any  other  programs. 

7.0  TEST  METHODS  AND  RESULTS 


The  following  program,  CHECK,  illustrates  the  use  of 
the  routines. 

CHECK  begins  by  setting  up  an  arbitrary  5x6  mesh  using 
a  data  statement  to  define  the  X  and  Y  arrays.  Next  a 
third  order  two  dimensional  polynomial,  11(1,  J) ,  having 
arbitrary  coefficients  is  evaluated  at  each  of  the  mesh 
points.  Since  NDFAULT  is  initially  set  to  zero,  CHECK 
is  required  to  supply  the  normal  derivatives  along  the 
boundaries.  The  equations  used  in  statement  300,  301  , 
and  302  were  obtained  by  differentiating  the  polynomial 
U (I , J) .  Next  subroutine  BICUB1  is  called  to  complete 
the  P,  Q,  and  S  arrays  and  the  results  are  printed. 
Following  this,  30  random  coordinates  are  generated 
over  the  x-y  mesh  using  a  uniform  random  number 
generator  available  in  the  CDC-3800  system  1 ibrary 
and  the  polynomial  U(I,J)  is  evaluated  at  each  point 
(IfFXACT) .  Also  subroutine  BICUB2  is  called  to  inter¬ 
polate  a  value  at  each  point.  Then  the  x-y  coordinates 
of  the  30  points  arc  printed  out.  Since  the  arbitrary 
polynomial  is  third  order  in  x  and  y,  the  interpolated 
and  exact  values  should  be  the  same.  The  fact  that  the 
elements  of  the  difference  column  are  essentially  zero 
(neglecting  rounding  errors)  indicates 
that  this  is  indeed  the  case.  Next  the  P,  Q,  and  S 
arrays  are  zeroed  and  the  parameter  NDFAULT  is  set  to  1, 
indicating  that  BICUB1  should  determine  the  edge 
conditions  rather  than  assume  they  are  supplied  by  CHECK, 
and  the  above  procedure  is  repeated.  Again  the  differ¬ 
ence  column  is  essentially  zero.  In  addition,  the  P,  Q, 
and  S  arrays  determined  by  BICUB1  agree  almost  exactly 
with  those  obtained  by  CHECK. 
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PR80RAM  CHICK  7170 

7180 

SAMPLE  MIOMN  T8  ILLUSTRATE  TH  USE  17  THE  7880 

BICUBIC  SPLINE  INTERPIlATJFN  PACNAOS  7|60 

7810 

SETTING  UP  «N  ARBITRARY  hESH  7920 

7990 

DATA<N«9>|<M*8),<MAXt6>  7940 

0ATA<X4t,,2,9,2,79,3,,9,),<Yn,il,9.2,2S',4,9,S,,7,3)  7990 

7980 

0 1 M| N S I |N  X(9),V<8), 1(9.1), P(9,6>, 0(8.9), $(8.9)  7970 

7980 

set  the  dimension  er  eaol  er  the  seven  mbrk  areas  to  max,  7990 

8000 

DIMENSION  Hl(6),W2(6),h3(6t,h4(«),u9(0)|Hl(0).W7(6)  8010 

8020 

evaluate  an  arbitrary  pblyninial  *r  Each  mesh  point  roso 

8040 

01  12  lil.N  8090 

Xl*x(|)  8080 

Dl  20  J«1«M  80*0 

V I  ■ V ( J )  8080 

U(I, J)a3,4l9,*Xl*17,«XI**2»l3,*Xt**3  8080 

X»Y!«(49?*28, tXI*ll,«Xl**2*l9i*XI**9)  8100 

X  *VIM2*(34,.6,»XI*1S,*XI»*!*43,»X!**3>  Bilo 

X  *Yla*3*(47,*2l,*Xl*l9,«X!**2*2i*Xl**3I  8j*0 

20  CONTINUE  8190 

12  CINTINUE  8140 

8190 
8180 

IF  NDFAULT  IS  SET  TO  0  ThE  ECQE  CONDITIONS  HILL  BE  INPUT  T§  0ICUI1,  8170 

IF  NDFAULT  IS  SET  TO  1  0ICLB1  h ILL  CALCULATE  THE  EDOE  CONDITIONS,  8180 

8180 

NCFAULTaO  8200 

8210 

CALCULATE  EXACT  EOQE  CENCITICNS  F|P  TEST  EQUATION  8220 

8290 

Nhl|N«l  8240 

M81|M«1  8290 

Dl  300  Ml.N.NMi  8280 

X|*X< I )  8270 

*2*0 

OET  N|8hAL  DERIVATIVES  HJTh  RESPECT  T|  X  8290 

8300 

Dl  300  J*1.M  8310 

YJ«y(J)  8320 

300  PH. J)*i9,*34,*Xl*249,»XI**2  8390 

X  ♦Yl«(j8,*38,«Xl*97,*Xl**7)  8340 

X  »VI**2*(8.*2I,*XI»129,*XI**2)  8390 

X  ♦Vl,,3*|2i,*30,*Xl*8,*XI*«2)  8380 

C  8370 

8380 

C  OET  NIRHAL  DERIVATIVES  H|Th  RESPECT  T|  Y  ALONG  EDGE  8990 

C  8460 

Dl  301  jal.M.MMl  8410 

Y I  ■ Y  f  J )  8420 
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Dfi  301  l*l»N  0430 

X|*XI1»  0440 

301  0  ( J  i !  )*40|*20i*Xl  ♦ie,»XlM2*i»,»X!Mj  0490 

X  «2l*ytM34i*6t*Xl*l3i«Xt«*2«43l*XfM3)  0400 

X  *3,*yl**2*(47i*21l*Xl*19l*Xt«*2*3l*Xl**3)  0470 

0400 

OET  NORMAL  DERIVATIVES  HJTH  RESPECT  T|  MOTH  X  AND  y  AT  EACH  CORNER  Of  04*0 
MESH,  0900 

0910 

DO  303  I«1«N,NM1  0920 

XI«X(J>  0930 

DO  302  J*liH|MMi  0940 

VI*Y< J)  0990 

302  3 < J « 1 l*20,*30l*Xl*97l*Xl*«3  0900 

X  *2i*VI*(A,*20i*X!*l20l*Xl**2)  0970 

X  #S|*yI*»2*(2li*30,*Xl*6,*XJ*«?)  0900 

09*0 

0000 

ESTIMATE  THE  AMOUNT  OF  T J Mi  RE6UIRID  FOR  A  CALL  TO  0ICUI1  0610 

0420 

200  TiTiMfiLEFTdl  0630 

0640 

COMPLETE  the  P.0,  and  S  ARRAYS  ,  THAT  is,  DETERMINE  NORMAL  DERIVATIVES  0690 
AT  EACH  MESn  POINT,  06A0 

0670 

•  . . . . . ••••*•*••**•••*•• . * . *****  0600 

CALL  BIcUBl(N,M,N0FALLT«Xiy«L|P,0«S(MAX|Wl«H2,H3,W4,H9iH4|U7|  ,  06*0 

•  *•*••••«••« . . . *•••  0760 

07lO 

0720 

CALCULATE  TIME  SPENT  |N  MILLISECONDS  0730 

0740 

T«(T*TImELEFT<1)>*1000,  0790 

PRINT  3  0760 

DC  100  1*1 1 N  0770 

DO  110  j*iiM  0700 

print  iiii!.jiX(n',Y(j)'tu,,n,R(3,,'»,Q(j,n,s<j',n  07*0 

110  CONTINUE  0000 

100  CONTINUE  0010 

PRINT  4*T  0020 

PRINT  2  0030 

0040 

use  system  random  number  qeneratcr  to  ooso 

GENERATE  RANDOM  COORDINATES  OVER  THE  X*T  PLANE  0060 

0070 

DO  10  K|l|30  0000 

Xl*RANF|«l>*(X(N J»X (1 ) )*X<i  )•  ,3*RANF<*1 J  00*0 

YL"iANrf«i)*4YtM»»VC1>>*YCi > • ,  3 * 3 ANr < «t >  0900 

UEXaCT*1,*19,*XI*17,*XI**2*03,*xI**3  0*10 

X*YI*(497*20|  *Xt*10  ,*XI**2*19»*XJ**3)  0*20 

X  «ViM2*(34,*«,*XI*13|»XI**2*43l*XI*»3)  0*30 

X  ♦Yl**3*|47,*2i,*Xl*i9,*XI**2*2«*Xl**3>  0*40 

0*90 

0*60 

ESTIMATE  THE  AMOUNT  OF  TIMS  REQUIRED  FOR  THE  CALL  To  BICU82  0*70 

0*00 
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FTN9.9A 


wimileftiii 

1NTERR0LATI  AT  VALUE  AT  (XIiYl) 


CALL  •lCUI2JXI,YI,LlNT,N|RR|R,N,M,)|fy(u,P,0,S) 


T»<?-T!HEL6fT(l)»*i0fl0, 

DIRf ■UInT*USXACT 

PRINT  1|K|X1 • V ! |UEXACT|UJNT,CIRF»N|NR|R|T 

10  continue 

see  if  he  have  completed  the  sf cend  pass 

IFINDFAuLT  ,E0,  1)  STIR 

REDCriNE  NDPAULT  80  THAT  THE  ECUE  C»ND!T!IN5  HILL  BE  CALCULATED 
AND  ZER|  OUT  THE  R«fii  ANC  8  ARRAYS, 

NDRAULTil 
DO  901  J*i,M 
DO  900  |*li N 

R(i,j)»o(j,n«s(j,i  ).o, 

900  CInTIhUe 

901  CONTINUE 

REPEAT  ASSUHlNQ  BOUNDARY  CONDITIONS  ARE  N|T  GIVEN, 

HERE  HE  USE  THE  LAOrANGE  INTERPOLATING  POLYNOMIAL  TO  DETERMINE  THE 
ED0E  CONDITIONS, 

00  TO  200 

R0RMAT  STATEMENTS  . •  «••* 

1  rORHATOlOX*  I9i2Xi4cri2l4»2»)*El6i9,«SXi  I9',2X,P18,4) 

2  rSRHATIiX^/, 

X  14X#*K*i9X,»XI*,l2X,*Yl*,0X»*UE*<ACT*|8X’,*UlNT*,l2X, 

X  •DirriR2NCE«.9X. TERROR*, SX<«TIHE*»/I 

3  format <ix(//,3xl*i*,3x#*j*,9xl*xM8xf*Y*,i2x,*u*,i2x,*P*,i2x,*Q*# 
X  l2X#*R*./> 

4  R0RmAT(i0Xi*THE  Call  TO  OICUBl  TOOK  APPROX, •,ri2,4,»  MSEC,*) 
in  R0Rm*TIiX,2< 13 ,1X 1 ,0(El2,4 , IX ) ) 

END 


9*0 


020 

0S0 

040 

090 

0*0 

070 

000 

0*0 

100 

110 

120 

130 

140 

190 

1*0 

170 

100 

1*0 

200 

210 

220 

230 

240 

290 

2*0 

270 

200 

2*0 

300 

310 

320 

330 

340 

390 

3*0 

370 

300 

3*0 

400 

410 

420 
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8..0  REMARKS 


Although  the  Lagrange  interpolation  procedure  described 
above  for  obtaining  the  "required"  boundaries  derivatives 
has  been  implemented  and  works,  a  different  and  possibly 
better  approach  [9]  to  this  problem  consists,  in  the 
analogous  one  dimensional  case,  of  not  having  a  break¬ 
point  [6]  at  the  second  and  second  to  last  data  points. 
Implementation  of  this  latter  approach  is  left  as  an 
exercise  for  the  interested  reader. 

For  those  readers  who  may  be  interested,  a  completely 
independently  conceived  and  different  set  of  ALG0L 
procedures  for  bicubic  spline  interpolation  is  given  in 
references  [10  and  11]  (in  German).  Comparison  of 
Spath's  algorithms  with  those  described  earlier  is  left 
as  another  exercise. 

As  a  further  remark,  it  should  be  obvious  that  the 
procedure  described  in  this  report  could  be  readily 
generalized  to  N  -  dimensional  cubic  spline  interpolation, 
where  N  is  greater  than  2. 
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